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Abstract 

An algorithm for parallel computation of transient response for structures is presented in which 
responses of substructures are computed independently for dozens of time steps at a time, and these 
substructure responses are then corrected to obtain the response of the overall coupled structure. 
The correction of the uncoupled substructure responses only requires the responses computed for 
interfaces at occasional points in time, and is done independently for different substructures in a 
very efficient procedure. A numerical example is presented to demonstrate the method and show the 
accuracy of the method. 

Introduction 

A significant amount of effort has been directed recently toward the development of methods for 
subdividing the computational effort associated with the solution of large transient response problems. 
The general approach of subdividing the computation associated with a given problem on the basis of 
a subdivision of the problem domain into subdomains has come to be known as domain decomposition 
in the last few years. 1,2 For transient response problems in structural dynamics, some efforts in this 
direction have been motivated by the need to solve problems for systems consisting of two or more 
well-defined subsystems, such as the Shuttle orbiter and its payloads, using modal data that have 
already been obtained for each of the subsystems rather than computing new modal data for the 
combined system. 3 ” 5 Other work has been done in the context of the element- by-element approach 
to finite element analysis. 6,7 More recently, Ortiz et al . have proposed methods specifically intended 
for concurrent computation of transient response based on a subdivision of the problem domain into 
subdomains. 8,9 In their approach, an implicit integration scheme is used to obtain response for each 
subdomain for a given time step, and the results of these computations are averaged at interfaces 
to yield an approximation of the response of the overall system. Hajjar and Abel have investigated 
the accuracy of these methods for certain structural dynamics transient response problems, and have 
concluded that their accuracy is inadequate for these problems when practical time step sizes are 
used. 10 

In all of the transient response methods mentioned above, computation of response on the sub- 
structure level can only be done independently for one time step at a time. In contrast to this, an 
algorithm was presented recently by these authors which allows independent computation of sub- 
structure response for an arbitrary number of time steps at a time. 11 After independent substructure 
responses have been computed, they are corrected based on the interface motion computed for sub- 
structures at each time step, to obtain the response of the combined structure. Allowing the response 
to be computed independently for a number of time steps at a time reduces the interdependence be- 
tween processors assigned to different substructures significantly, which can be important when the 
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amount of computation required for different substructures is unequal. Also, if there are more sub- 
structures than processors, the cost of swapping different substructures in and out of processors will 
be reduced if it can be done less frequently. 

In the present paper, an extension of the algorithm presented in Ref. 11 is presented in which 
independent substructure response computation can proceed for much longer periods of time. Inde- 
pendent substructure responses are corrected on the basis of computed interface motion sampled at 
occasional points in time. The correction procedure for obtaining the response of the structure from 
the computed substructure responses is extremely efficient once the transient response computation 
is under way, although there is some computational overhead required to set up the correction capa- 
bility. A numerical example is presented which illustrates the method and shows the accuracy that 
is obtained. 


A Method Using Substructure-Level Response Computation 

The algorithm presented in this paper is for computing the transient response of structures whose 
motion is governed by the equation 

Mu + Cii + Ku = F(t) (1) 

where Af, C, and K are taken to be constant mass, damping, and stiffness matrices, ii, u, and u are 
acceleration, velocity, and displacement vectors, and F{t) is a vector of forces exciting the system. 
As mentioned in the Introduction, the transient response of a given structure is computed in this 
algorithm by solving transient response problems for the substructures defined by decomposing_the 
structure. To introduce the notation that will be used in this paper, a mass matrix for a structure 
composed of two substructures is shown below, after a possible reordering of rows and columns: 


M = 




0 

*4! 

M.S + Ml 5 

"a 

0 

Mls 

<1 


( 2 ) 


The superscripts in parentheses tell which substructure a given matrix partition is associated with, 
and the subscripts S and L refer to matrix partitions associated with shared , or interface, and locals 
or internal degrees of freedom. For some of the development in this paper, a structure composed of 
only two substructures is considered in an effort to simplify the presentation. However, the methods 
presented will be applicable for an arbitrary number of substructures. 

Because responses will be obtained for each of the substructures a structure is composed of, some 
convention must be adopted for representing the structure response in terms of the substructure 
responses, particularly at the interfaces. In this paper, the approach taken is similar to the standard 
approach for the assembly of element matrices in the finite element method. The response of the 
structure in interface degrees of freedom is represented as the sum of the interface responses for the 
substructures sharing the interface, e.g., 


u = 
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( 1 ) 
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so that each substructure’s interface response is only one component of the total interface response 
of the structure. Of course, if this convention is adopted, substructure transient response problems 
must be defined and solved in such a way that the response of the structure obtained by assembling 
together the substructure responses is accurate. 
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Substructure response problems can be defined for independent computation by extracting equa- 
tions from the structure equations of motion, and they will be of the form 



where “hat” symbols identify matrix or vector partitions for which a policy for assigning the cor- 
responding partitions in the structure equations of motion to the different substructures must be 
determined. Again, reordering of rows and columns may be necessary to collect all “shared degrees 
of freedom together for a given substructure. Simply computing substructure responses that satisfy 
these equations and assembling them together will not result in an accurate representation of the 
response of the overall structure, because the interaction between substructures is neglected in such 
an approach. It must be noted that in the response of the structure, each substructure has two 
sources of excitation. One is the external applied force, which appears on the right hand side of the 
equation above, and the other is due to interaction with adjacent substructures at the interfaces. 
This suggests a two-step approach for computing the responses of substructures in the response of 
the coupled structure. The first step consists of obtaining independent substructure responses that 
satisfy the substructure equations of motion above. These responses neglect any interaction between 
substructures. Then the second step consists of correcting these substructure responses to obtain 
responses of substructures in the motion of the coupled structure. It will be shown that this second 
step can be accomplished with a surprisingly small amount of effort, and with very little information 
from the independent substructure responses, 

If independent responses satisfying the substructure equations of motion are computed, and as- 
sembled together and inserted into the structure equations of motion, a residual r(t ) will be obtained. 
For a two-substructure structure the residual will be given by 
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By making use of Eq. (4), the residual can be obtained as 
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+ (4s + 4 2 i)(4 1) + 4 2) ) - 4s4 1) - 4 2 J4 2) 

+ (43 + 4 2 i)(4 1) + 4 2) ) - 4 1 J4 1) - 4s4 2) 

-Fs + Ff + F™. ( 7 ) 


Note that the residual associated with one substructure is given entirely in terms of the interface 
motion computed for adjacent substructures. Note also that rs(f) is defined in terms of the “hat” 
partitions of Eq. (4), and can be obtained as a null vector, if these “hat” partitions are chosen to 
satisfy the following: 
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With this as motivation, the “hat” partitions are taken to be defined this way in this paper. A physical 
interpretation of this choice is that for each of the independent substructure response problems, the 
structure is modeled as if it were clamped one node beyond the interfaces, and the excitation acting 
on the structure at the interfaces is divided between the substructures that share the interfaces. 

The residual in the equations associated with a given substructure can be seen to be a result 
of including the interface motion of adjacent substructures in the given substructure’s equations 
of motion. This interface motion for adjacent substructures was neglected in the solution of the 
independent substructure response problems. In order to obtain the true response of the structure, 
the substructure responses must be corrected to account for adjacent substructures’ interface motion, 
so that when the substructure responses are assembled into the structure equations of motion, the 
residual is zero. 

For the correction to the first substructure’s response, note that if the interface motion for the 
second substructure were given, the residual in the structure equations of motion associated with 


the first substructure would be defined. The first substructure’s response would have to be corrected 
by adding a response of the first substructure to the negative of the residual resulting from the 
interface motion of the second substructure. The second substructure’s response would have to be 
corrected in a similar manner, if the interface motion for the first substructure were given. However, 
the interface motion for both substructures is not known a priori , because all of the interface motion 
will be changed as a result of the corrections to the substructure responses. The responses of both 
substructures will have to be corrected simultaneously, so that the response of each substructure to 
the negative residual due to the other’s corrected interface motion will be addecTto the independently 
computed substructure response. The following paragraphs present a method for accomplishing this. 

Because the residual is defined in terms of interface motion, it is convenient to introduce a vector 
containing the interface accelerations, velocities, and displacements for the &th substructure 
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With this definition, the correction of the first substructure’s response to account for the second 
substructure’s interface motion will be the response to an excitation of the form 


f m (> ) = 


-Mg 




( 10 ) 
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where the degrees of freedom are ordered as in the structure equations of motion. If the second 
substructure’s interface motion is given only at the beginning and the end of a time interval 

consisting of p time steps of length At, the interface displacement can be approximated over 

the time interval 0 < t < pAt by interpolation. Hence, u^\t) is assumed to take the form 


4 2) (0 = MW MM]^ 2 \ 0) 

+ imm mm mmM 2) {pa t) ( ii ) 


where I represents a unit matrix and i>i(t), i = 1, ... ,6 are interpolation functions that must satisfy 
the following end conditions: 


MO) = i, 
V>2(0) = i, 
MO) = i, 

MpAt) = 1, 
MpAt) = i, 
Mp At) = 1 , 
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Mp) = M(o) = MpAt) = Mp^t) = MpAt) = o, 
^4(0) = ^4(0) = ^4(0) = MpAt) = MpAt) = 0, 
Mo) = M(0) = ^ 5 ( 0 ) = MpAt) = MpAt) = 0 , 
M(0) = ^ 6 (0) = ^6(0) = MpAt) = MpAt) = 0. 
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Quintic polynomials were used for the results obtained in this paper. Expressions for ii s (t) and 
iSg\t) for defining the excitation for correcting the first substructure’s response are easily obtained 
by differentiating the interpolation functions. 

With vfg\t) defined in terms of u< 2 )(0) and vW(pAt), the corrected interface motion for the 
first substructure at the end of the time interval will be the sum of the response to the independent 
response problem and the response based on r>( 2 )(t), 0 < t < pAt. Hence, it will have the form 


v^(pAt) = u^j(pAt) + Si 2 t>^(0) + T\2v( 2 \pAt), (13) 


where each column of the matrices S 12 and T 12 contains the first substructure’s interface response 
at t = pAt to a negative residual specified by a column of the first or second matrix, respectively, 
on the right-hand side of Eq. (11). Using a similar approach, the corrected interface response of the 
second substructure at the time t = pAt can be expressed in terms of the first substructure’s interface 
motion as 

t>( 2 )(pA<) = t>£jr(pA*) + S 2 it> (1) (0) + T 2 M l \pAt). (14) 


As mentioned above, corrected interface' motion for an adjacent substructure is not known before 
the reconciliation is accomplished. All that is known in the two equations above is the interface 
motion of both substructures at t = 0, from initial conditions, and the interface motion obtained 
from the solution of the independent substructure transient response problems. However, given the 
set of linear equations in Eqs. (13) and (14), it is straightforward to solve for the unknowns, with the 
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More compactly, the reconciled interface motion is given by 
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where the matrices S and T are readily identified. The corrected motion for the first substructure’s 
local degrees of freedom at t = pAt is given by 


/ u L\p At ) 
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(17) 


where columns of the matrices 5^ and T 'jp contain responses in local degrees of freedom to interpo- 
lation functions for representing interface motion. These two matrices are naturally obtained at the 
same time that the matrices 5i2 and T 12 are obtained, from the solution of the same substructure 


response problems. The corrected motion in local degrees of freedom for the second substructure is 
obtained in the same manner. Once the motion in both local and shared degrees of freedom has been 
corrected for t = pAt , the initial conditions have been obtained for ongoing computation of response 
for the next p time steps. 

The developments presented here are easily applied to structures composed of more than two 
substructures. For example, if there are three substructures, the matrices S and T in Eq. (16) take 
the form 
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and modification of the rest of the procedure presented for two substructures is straightforward. 


Infrequent Reconciliation of Substructure Responses 

In the method of the preceding section, responses are computed independently for different sub- 
structures for p time steps at a time, and then the independent substructure responses are corrected 
to obtain substructure responses in the response of the overall coupled structure. In this section, a 
procedure for carrying out the reconciliation of independent substructure responses after a number of 
p-step time intervals is developed. This procedure will allow substructure responses to be computed 
independently for long periods of time without correcting for interaction between substructures. 

The interface motion for the second substructure over the time interval pAt < t < 2 pAt can be 
approximated in terms of the interpolation functions introduced in the preceding section and the 
interface motion at the beginning and end of the time interval as 

4 2) (0 = [Mt*)I V»3(C)/]u< 2 )(pA0 + [^(i*)/ j 5 (t*)I Mnnv {2) (2pAt), (19) 

where t* = t — pAt. Recalling that substructure responses have two components including the 
response to external excitation, which is represented in the independent substructure responses, and 
the response due to interaction with adjacent substructures, which is represented in the correction to 
the independent substructure responses, the interface response of the first substructure at the time 
t = 2pAt will have the form 

v^(2pAt) = 4£(2pA<) + 5i 2 (2pAO® (2) (0) + T u (2pAt)v^(pAt) + T 12 (pAt)v^(2pAt). (20) 

Here, the columns of SuftpAt) contain responses of the first substructure at t = 2pAt based on 
the second substructure’s interface motion, which is given in terms of the interpolation functions ip\ , 
ip2, and ip 3 for 0 < t < pAt, and is extended as zero for pAt < t < 2pAt. Similarly, the columns 
of Tu(2pAt) contain responses of the first substructure at t = 2pAt based on interface motion of 
the second substructure which is given in terms of the interpolation functions V’s? and for 
0 < t < pAt, and is extended in terms of ip i, ip 2 , and ip 3 for pAt < t < 2pAt. The matrix Ti 2 (pAt) 
is simply the matrix Ti 2 of the preceding section. 
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The interface motion for both substructures at t = 2pAt can be written as 
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with 5" and T matrices defined in terms of 0, S\ 2 , S 2 i, 0, etc., as in the last section. Solving for 
®W(2pAi) and t?^(2pAf) gives the result 
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and letting A = [7 - r(pA<)] _1 , 5, = S(ipAt), and T; = T(ipAt), v(2pAt) can be obtained in terms 
of initial conditions and independent substructure responses as 

f ®(°) } 

v(2 P At) = A[(S 2 + T 2 AS!) T 2 A 7] < v ind (pAt) V. 

I v ind (2pAt) J 

The corrected interface motion at t = 3pAt can be found using the same approach. When the 
interface motion for the different substructures is assumed in terms of interpolation functions as in 
Eq. (19), linear equations involving u(3pAi) can be written as in Eq. (21). These equations can be 
solved for v(3pAt), yielding the result 

r ®(°) 

„(3pAl ) = M[S 3 Tj T, i J] ■ (26) 

l Vi nd (3pAt ) , 

Interpolation functions are simply extended as zero into the time interval 2pA< < t < 3pAt in the 
generation of responses for matrices S3 and T3. Inserting the expressions for v(pAt) and v(2pAt) 
from Eqs. (24) and (25) gives v(3pAt) in terms of initial conditions and independent substructure 
responses as 

v(3pAt) = A[(S 3 + T 2 AS 2 + (T 3 A + (T 2 A) 2 )S 1 ) (T 3 A + (T 2 A) 2 ) 

®(0) ] 

Vind(pAt) I 

fmd(2pAt) | ‘ 

y Vi nd (3pAt) J 

This result can be generalized for finding the corrected interface motion at a time t = mpAt, with 
the result that 
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where Bq = J, and the other Bi matrices are defined by the recursive formula 


t-i 


Bi = J2(T,- l+ iA)B h 

1=0 


(29) 


so that Bt = T 2 A, B 2 = T 3 A + ( T 2 A ) 2 , B 3 = T 4 A + T 3 AT 2 A + T 2 AT 3 A + ( T 2 A ) 3 , etc. Defining a 
matrix C m as 

Cm = a \ ( B m - 1 B m - 2 *'• , (30) 

^ »=0 / 


the corrected interface motion can be obtained separately for each substructure by partitioning C m 
into upper and lower halves Cm^ and Cm \ and multiplying each by the vector on the right hand side 
of Eq. (28)* For parallel computation, if different processors are assigned to different substructures, 
the processor for the kth substructure only needs to have access to Cm ^ and the interface motion 
computed independently for all substructures for every pth time step. 

After interface motion has been corrected for t = mpAt, the motion for local degrees of freedom 
for each substructure can be corrected. As an example, the corrected local motion for the first 
substructure will be given by 


f u^\mpAt) If u { l] nd (m P At ) 
1 ii ( l\mpAt) j \ u ( H d (mpAt) 
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where the matrices S ^ and contain responses in local degrees of freedom to interface motion 
given in terms of interpolation functions, and are analogous to the Si and T t matrices used above in 
terms of subscript numbering. The vector of the second substructure’s corrected interface motion at 
every pth time step is given in terms of the independently computed interface responses as 
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Therefore, the product of the matrix on the right hand side of Eq. (31) and the matrix on the right 
hand side of Eq. (32) is the matrix by which the vector of independently computed interface responses 
must be multiplied to obtain the correction for the motion in local degrees of freedom for the first 
substructure. The same approach is taken to find the correction for the motion in local degrees of 
freedom for the second s ubst ruct ure. 

To summarize, the developments presented in this section permit the independent computation 
of response for different substructures for a total time interval of length mpAt. The interface motion 
for all of the substructures at the end of this time interval can be corrected using Eq. (28), and then 
the motion for local degrees of freedom for each of the substructures can be corrected as shown above. 
Once these corrections are made, initial conditions are obtained so that independent computation of 
substructure responses can proceed again for another mpAt. The amount of computation required 
for the corrections is very small compared to the amount of computation required for obtaining the 
independent substructure responses. The computational “overhead” that is required for this method 
consists of obtaining substructure responses to interface motion specified in terms of interpolation 
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Figure 1: Plane truss used in the numerical example, and its division into substructures. 

functions, and carrying out the matrix operations outlined above to obtain the matrices required for 
making corrections. This overhead is justified if the transient response of the structure must be com- 
puted for a long time. The amount of computation required both for the “overhead” operations and 
for the corrections is determined by the dimensions of the matrices involved, which is determined in 
turn by how many shared and local degrees of freedom are associated with each of the substructures. 

Numerical Example 

The algorithm presented in this paper is demonstrated on an example structure which is shown 
in Fig. 1. The structure is a plane truss composed of 143 aluminum members, each of which has an 
elastic modulus of E = 70 x 10 9 N/m 2 , a cross-sectional area of A = 4 x 10“ 4 m 2 , and a density of 
p = 2710 kg/m 3 . The dimensions are as shown. A force is applied to the top right corner of the truss 

starting at t = 0, and it is given by 

F(t) = 5(1 — cos f It) (Newtons), (33) 

where fl = 590.3 radians per second, which is between the second and third natural frequencies of 
the structure. The truss has eighty-eight degrees of freedom, and is assumed to have proportional 


533 



6 



xlO- 6 



0 0.02 0.04 0.06 0.08 “ 0.1 0.12 


Time t, sec. 


Figure 2: Plots of exact response (dashed line) and computed response (asterisks). 


damping of the form C — otM + f3I { , where a and /? are chosen to give modal damping factors 
between one and five percent. For application of the algorithm presented in this paper, the structure 
was partitioned at the top of the sixth bay into two substructures, which are also shown in Fig. 1. 
Note that each substructure is modeled in the algorithm as being effectively clamped one truss bay 
beyond the interface, as shown in the figure. 

In Fig. 2, the horizontal displacement of the structure at the point where the excitation is applied 
is plotted. The dashed line is a plot of the exact response, obtained from a mode-by-mode exact 
solution, and the asterisks represent values that were obtained using the algorithm of this paper. The 
responses of the two substructures were obtained using an algorithm that finds the exact response to a 
piecewise linear approximation of the excitation. 12 A time step of At = 3.74 x 10 -4 seconds was used, 
which is equal to about one twenty-eighth of the period of the excitation, and is also approximately 
equal to the period of the highest mode of the structure. For larger time steps, the error becomes 
visible on a plot scaled as in Fig. 2, when the piecewise linear algorithm is used on the structure 
as a whole. In this example, substructure responses were computed independently for sixty time 
steps at a time, and then corrections to the independent substructure responses were made based 
on the interface motion computed for every tenth time step. Therefore, the quintic interpolation 
polynomials for interface motion were defined over time intervals of length pAt with p equal to ten, 
and there were six of these time intervals in each time period over which independent substructure 
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responses were computed. 

Because the response of the structure was only corrected for every sixtieth time step, the asterisks 
on the plot in Fig. 2 are sixty time steps apart. It should be noted, however, that the response for 
any degree of freedom at any time can be obtained in a straightforward manner with a small amount 
of additional computation. From the plot of Fig. 2, it is evident that the accuracy obtained in this 
example is quite adequate for most purposes, even though the corrections to independent substructure 
responses were made based on a very limited amount of information. The only approximations 
made in obtaining these results were in the piecewise linear approximation of the excitation and the 
piecewise quintic approximations of the interface motion. 

Summary 

In this paper, an algorithm is presented for computing the transient response of structures by 
computing the transient responses of substructures. The algorithm is well suited for parallel im- 
plementation, where a different processor would be assigned to each substructure. The fact that 
computation can proceed independently for different substructures for dozens of time steps at a time 
reduces the interdependence between processors, which can be of considerable importance when dif- 
ferent substructures require different amounts of computational effort per time step. The correction 
of independently computed substructure responses to obtain the response of the structure acting as 
a whole requires only the interface jnotion computed for substructures at occasional points in time. 
This correction of substructure responses can be done independently for different substructures once 
the interface motion for all of the substructures has been computed, and this correction requires 
very little effort. Because of this, the total amount of computation required using this approach will 
be only slightly greater than the amount required to solve the transient response problem for the 
structure as a whole for many problems. A surprisingly high level of accuracy is obtained using this 
algorithm, in view of how little information is required for making corrections to the independent 
substructure responses. 
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